Reconstruction of spatially inhomogeneous dielectric tensors via optical tomography 
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A method to reconstruct weakly anisotropic inhomogeneous dielectric tensors inside a transpar- 
ent medium is proposed. The mathematical theory of Integral Geometry is cast into a workable 
^+ ■ framework which allows the full determination of dielectric tensor fields by scalar Radon inver- 

C^ ' sions of the polarization transformation data obtained from six planar tomographic scanning cycles. 

C^ , Furthermore, a careful derivation of the usual equations of integrated photoelasticity in terms of 

^S) . heuristic length scales of the material inhomogeneity and anisotropy is provided, making the paper 

a self-contained account about the reconstruction of arbitrary three-dimensional, weakly anisotropic 
dielectric tensor fields. 
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I. INTRODUCTION 



The inverse boundary value problem of recovery of anisotropic spatially varying electromagnetic properties of a 

O , medium from external measurements at a fixed frequency is among the most mathematically challenging of inverse 

O ■ problems. For low-frequency electromagnetic measurements where a static approximation is valid, an anisotropic 

^ , dielectric tensor or conductivity tensor is uniquely determined by complete surface electrical measurements up to 

O ' a gauge condition [1]. For isotropic [2] and chiral isotropic media [3] a knowledge of all pairs of tangential electric 

' c/) ' and magnetic fields at the boundary, for a single non-exceptional frequency, is known to be sufficient to recover the 

^~>. material properties. Anisotropic materials are important in many practical problems with anisotropy arising from, 

T^ ' for example, flow, deformation, crystal or liquid crystal structure, and effective anisotropic properties arising from 

, M^ , homogenization of fibrous or layered composite materials. For general anisotropic media the question of sufficiency of 

data for reconstruction remains open. 
fS| , In this paper we concentrate on a specific high-frequency case of considerable practical importance. We assume that 

t> ■ the material is non-magnetic, i.e. has a homogeneous and isotropic permeability equal to the vacuum value ^o; that 
^^ \ the conductivity is negligible; and that the permittivity, or dielectric tensor, is weakly anisotropic in a sense that will be 
defined below. We present a mathematical framework which allows us to invert the polarization transformation data 
obtained from tomographic measurements of light rays passing through an optically anisotropic material at sufficiently 
(^ , many angles for the six independent components of the (spatially varying) dielectric tensor gy inside the specimen. It 
^+ ' will be shown that, for weak anisotropy, our method allows the full determination of all tensor components, provided 
f^ , that tomographic measurements are made for six carefully chosen orientations of the planes in which the light rays 
^^ ' scan the medium. 

O \ The equations describing the passage of light through inhomogeneous and weakly anisotropic media have been for- 

' ^ ■ mulated in Refs. 4-6, and in section II we shall show how to obtain these equations from a geometric-optical starting 
^~j' point by expanding the electric field in powers of appropriate length scales which describe the inhomogeneity and 
T^ , anisotropy of the medium. Also, a condition for "weak" anisotropy will be specified upon which we shall linearize 
^H' (section IV) the equations determining the polarization transfer matrix along the light rays. In section V we then 
^ , present our method of reconstructing the permittivity tensor in the linearized limit by performing scalar Radon trans- 
forms on the polarization transformation data obtained from six scanning cycles on different planes intersecting the 
. , specimen. In section VI the method is demonstrated by giving a visual example of reconstruction of the permittivity 
j3 ■ tensor inside an axially loaded cylindrical bar: in this case, the stress tensor in the (transparent) medium gives rise to 
optical anisotropy, and it will be seen that our method works well in the limit of weak anisotropy. The mathematical 
basis for this method has been anticipated in the seminal work of Sharafutdinov[7], whose book contains a general 
theory of "ray transforms" of tensor fields on n-dimensional Euclidean spaces, and an examination into the possibility 
to invert them for reconstructing the underlying tensor fields. Our work presented here is an adaption, and par- 
tial reformulation, of this highly mathematical framework, to the specific objective of reconstructing inhomogeneous 
dielectric tensor fields via optical tomography. 

It should be mentioned that the determination of anisotropic dielectric tensors as presented here has a closely related 
variant in the problem of reconstruction of a stress tensor field inside a loaded transparent material; the phenomenon 
that an initially isotropic medium becomes optically anisotropic when under strain is called Photoelasticity [8-12] and 
may be used to obtain information about the internal stress in a strained medium from polarization transformation 



data obtained by tomographical methods. The photoelastic effect as a means to retrieve information about internal 
stresses has been studied extensively: A method termed "Integrated Photoelasticity" [13-17] is well-known, and it 
was pointed out[18-21] that information about the difference of principal stress components could be retrieved from 
appropriate Radon transformations of polarization transformation data. However, these methods do not succeed in 
reconstructing the stress components separately and therefore the full stress tensor, and the linearly related dielectric 
tensor, can be obtained in this way only for systems exhibiting a certain degree of symmetry, such as an axial 
symmetry. Other methods of reconstruction have been suggested, for example, a three-beam measurement scheme [22], 
where for axisymmetric systems an onion-peeling reconstruction algorithm was proposed[23, 24]. Another method, 
which in principle is capable of determining a general three-dimensional permittivity tensor, is the "load incremental 
approach" [25]. Here, the stress on the object is increased in small increments, and at each step, a measurement cycle 
is performed. 

The new results in this paper arc twofold: Firstly, we derive the set of equations determining the polarization transfer 
matrices for a given light ray scanning the object. These equations are somewhat related to the standard equations 
of integrated photoelasticity [14], but here we present a more rigorous exposition of the heuristic length scales [6], and 
the approximations based thereupon, which underlie the usual derivation of these equations from Maxwell's equations 
in a material medium; this is done in sections II, III and IV. Secondly, we present a novel scheme for reconstruction 
of arbitrarily inhomogcncous dielectric tensors in the interior of the specimen, subject only to the condition that the 
birefringence is "weak" in a sense which will be specified in sections V and VI. 

II. HEURISTIC LENGTH SCALES 

It is assumed that the material is non-absorbing for optical wavelengths, and furthermore is non-magnetic, fiij = 
fJ'oSij, where ^o is the magnetic permeability of vacuum. As for the permittivity we assume that the dielectric tensor 
deviates "weakly" from a global spatial average 

where B denotes the body having a volume volB, tr e — X]i=i ^a is the trace operation, and e — e(x) is the dielectric 
tensor. In a zero-order approximation the material therefore may be regarded as homogeneously isotropic, with scalar 
dielectric constant e as defined in (1). Typically, this behaviour of weak deviation from a homogeneously isotropic 
background will be satisfied for glasses and certain resins under moderate internal stress or external load. The scalar 
constant of permittivity e will be a reference quantity when we specify the dimcnsionless degree of anisotropy in 
eq. (8). 

For the actual problems relevant to our work, the length scales characterising inhomogeneities in the material 
are much larger than the wavelength of the (monochromatic) light passing through the object, so that the usage of 
geometric-optical approximations is justified. Hcuristically, two such length scales can be conceptualized [6]: A scale 
Iq characterising the degree of inhomogcneity may be introduced by 

;o|K-Vtre(x)|-|tre(x)| , (2) 

where k, denotes a unit vector in the direction of wave propagation. Furthermore, in any anisotropic medium, two 
preferred polarization directions ep(x), p = 1,2, for each given direction of wave propagation[26-30], exist at each 
point, so that a scale Ip measuring the rate of variation in these polarization directions can be speficied by 

^pIk- Vep(x)| -- |ep(x)| . (3) 

These scales should be compared to the average wavelength A of the monochromatic light passing through the medium, 
so that, when the limit 

Ip, ?o > A (4) 

is satisfied, the (complex) electric field may be given in the form 

E(x,t) =E(x)e^'^(^)-*'^* , (5) 

where the eikonal (f>{x + As) = 0(x) -I- As • V0(x) -I- O (A/Zq) describes a locally-plane wave with wave vector \/(f>{x), 
and the amplitude A k ■ VE ~ O (X/lp) E varies weakly on the length scale A. Motivated by these considerations. 



Fuki, Kravtsov, Naida (FKN)[6] introduce a dimcnsionless scale 

(^ ^] 
a = max i -Tt -r f ■ (6) 

The limit of geometrical optics then can be specified by the condition 

a < f . (7) 

We also need an explicit measure of anisotropy 

(3 ^ ma.x\\Aij\\ , (8b) 

where (3 is an appropriate number characterising the magnitude of the components of the dimensionless anisotropy 
tensor Aij, such as a global maximum, etc. If anisotropy is not weak, then at each point in the medium we have a 
continuous splitting of rays due to the fact that there are two distinct phase velocities, and two distinct ray velocities, 
for any given propagation direction. A condition for weak anisotropy therefore arises if we demand that the passage 
of light through the material can be described by a single ray which is influenced by the local variations of the 
optical tensors only in the way that the polarization directions rotate. This is the situation that commonly occurs in 
photoelasticity and is also of greatest interest to our work. 
It was shown by FKN[6] that ray splitting can be ignored, if 

-<1 • (9) 

a 

In this case [6] "it is impossible to discriminate experimentally between split rays", and one can effectively replace the 
two rays by a single isotropic ray which is obtained from the isotropic part of the dielectric tensor alone. This is indeed 
the domain we are most concerned with, since experimentally no ray splitting is seen in photoelastic experiments. In 
fact, for those applications we are interested in it is usually true that light propagates along straight lines through 
the specimen, so that the trial solution (5) may be replaced by an even stronger ansatz 

£:(x,t) = E(x)e*'^-''-*"'* (10) 

with constant wave vector k, just as for a plane wave. The phase velocity associated with this wave vector is the one 
associated with the average permittivity e defined in cq. (1), 

k^- , u = ^ , A = ^ . (11) 

However, the amplitude E(x) has a spatial dependence which accounts for the variation of the two preferred polar- 
ization directions along the light ray. 

Under the conditions (7) and (9), the electric field E and electric displacement D behave like 

K-n^of^jn , K-E-O^^+Z?) E . (12) 

This means that D is nearly transverse, while the same is true for E only if we assume in addition that the degree of 
anisotropy is small, 

/3<1 . (13) 

This is the condition of weak anisotropy, and our method is formulated for this "quasi-isotropic" regime. 

III. EQUATIONS SATISFIED BY THE TRANSFER MATRICES 

The information about the change in the state of polarization of a light beam passing through the material is 
encoded in a two-dimensional unitary transfer matrix. The equation satisfied by the transfer matrix along a light ray 
is given in most accounts on photoelastic tomography[18-24], but the various approximations taken in the process 



of neglecting higher powers of ratios A/Zq and \/lp are not always stated clearly; we therefore briefly summarize the 
necessary steps here: 

On inserting (10) into Maxwell's equations we obtain 



kx 
It is easy to show that 



kx E 



- u) /ioeE — i 



V X (k X E) + k X (V X E) 



- V X (V X E) = 



(14) 



V X V X E-O 



E 



(15) 



hence the term can be neglected in the geometrical-optical limit (7). Then (14) takes the form 



K;x(KxE)--^|vx(AtxE)+K;x(VxE)| + iiqu^eE = , 



(16) 



where u was given in eq. (11). The longitudinal component of eq. (16), obtained by projection onto the unit vector 
K, = k/A: in the direction of propagation of the light beam, is of the order 



Ol^lE 



(17) 



and is neglected in the geometrical-optical limit. Hence we only retain the transverse components of E and D, i.e. 
the components perpendicular to the wave propagation n. 

Now let us study eq. (16) in a coordinate system in which the direction of propagation k is along the z axis. 
Then (16) becomes 
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where Aij was defined in eq. (8a). The solution of (18) can be expressed via a transfer matrix 



Ei{z) 
E2{z) 



U{z,zo) 



Ei{zo) 
E2{zo) 



(18) 



(19) 



where U satisfies an ordinary differential equation similar to (18), together with initial condition 

d TT 

— U{z,zo)^i-Aj_{z)U{z,zo) , 
clz A 



(20) 



U{zq,zo) = I2 



1 
1 



Eqs. (20) can be expressed as an integral equation 



U{z, zo) = l2 + i- I dzi A^{zi) C/(zi, zo) 



(21) 



where A±_ denotes the matrix of transverse components of Aij as they appear in eq. (18). A formal solution of (21) 
is given by the Born-Neumann series 



C/(z,zo) == I2 + {i^ / dzi Ax{zi) + {i^ / dzi A^{zi) / dz2 A^{z2) + 



(22) 



The transfer matrix U is unitary and thus preserves the norm of the complex electric field vector. Physically this 
means that intensity is preserved, so unitarity here just expresses energy conservation of the light ray. This must 
indeed be the case, since we have assumed a non-absorbing medium. 



IV. THE LINEARIZED INVERSE PROBLEM 

Assuming that the transfer niatriees U have been determined for sufficiently many hght rays scanning the medium, 
the associated inverse problem now consists in reconstructing the anisotropy tensor Aij from the collection of these 
transfer matrices; this inverse problem is obviously non-linear in A^j, as can be seen from eqs. (21 and 22). The 
solution to the fully non-linear problem is not known as yet. However, in the quasi-isotropic regime we can deal with 
the linearized inverse problem: this is defined by a truncation of the Born- Neumann series in (22) after the first-order 
term 

U{z, zo) = ^2 + ij I dzi A^{zi) . (23) 

For example, for a relative anisotropy oi ly ^ 10~^, a wavelength of A ^ 0.5 x 10^® m and assuming a length of 
L ~ 1 m of the object, we find that the first-order term in (23) is of the order 10^'^, hence the linearization will be a 
good approximation in this case. 

The transfer matrices U must be determined by suitable measurements of the change of polarization along each light 
ray passing through the medium at many different angles, possibly including interferometric methods. By measuring 
three so-called characteristic parameters [13] the 5'C/(2)-factor of the transfer matrices U can be computed using 
the Poincare equivalence theorem[31], a matrix decomposition theorem which allows to interprete the characteristic 
parameters as the retardation angle A, the angle of the fast axis 9, and the rotation angle 6 of an equivalent optical 
model consisting of a linear retarder and a rotator with these optical parameters. The Poincare equivalence theorem 
can be formulated in terms of Jones matrices or Stokes parameters on the Poincare sphere; a contemporary exposition 
of these relations was given recently in Ref. 32. 

Thus, the measurement of A, 6, S determines a unimodular matrix 5'(A, 9, S) such that 

U^S{A,0,S)xe''^ , S{A,e,S)eSU{2) , (24) 

where $ is the global phase of the transfer matrix U . In the general case, assuming no restrictions on the degree of 
anisotropy, this global phase can be arbitrarily large, and furthermore cannot be determined from the characteristic 
parameters alone, but must be measured e.g. through interferometric methods, for each given light ray. However, 
in the limit of weak anisotropy, the global phase is effectively determined by the characteristic parameters alone: 
Suppose that the unimodular matrix S has been computed by measuring the parameters A, 9, 6 using the Poincare 
equivalence theorem; from eq. (23) it then follows that the unknown phase <& must be chosen so as to make the real 
parts of the diagonal (non-diagonal) matrix elements equal to one (zero), 

} ( } { (25) 



We can therefore determine the phase $ from any of these equations, or, for the sake of numerical stability, from all 
of them, so as to obtain an average value of $. Thus, in the weak-anisotropy limit, the number of real "degrees of 
freedom" in the transfer matrices is effectively equal to three, rather than four as in the general case. 

V. SOLUTION OF THE LINEARIZED INVERSE PROBLEM BY SIX SCALAR RADON INVERSIONS 

We now show that the linearized inverse problem can be reduced to six scalar Radon inversions performed on the 
polarization transformation data. We first specify a plane P(y, rj) in R^ which intersects the specimen and contains 
the point y; the orientation of the plane is determined by a unit vector t] normal to the plane. Consider the straight 
line t i-^ y + tK with unit vector k lying in P{y,ri), describing a light ray passing through the specimen and lying 
in the given plane -P(y, rj). Let ^ be a third unit vector, orthogonal to k and r) in such a way that (^, rj, k) form a 
right-handed system. Then the equation describing the polarization transformation along the light ray in the direction 
K is given by the analogue of (23), 






1 
1 



to 



A^^{ti) A,^r,{tl) 



(26) 



The measurement of characteristic parameters determines the transfer matrix on the left-hand side of (26). We now 
repeat these measurements for all lines lying in the given plane P(y, J7) and thus obtain a collection of line integrals 
for the normal component A^^ in P, 

dti A^^(ti) = \u^^{+oo, -00) - ll — (27) 



ITT 



for any pair of directions k and ^; we could extend the limits of integration in (27) to ±00, since in practical situations 
the object will be placed inside a tank with a phase-matching fluid, hence the value of Aij outside the object is zero. 
The set of line integrals in (27), taken for all light rays in P, is called the transverse ray transJorm\7] of Aij with 
respect to r). 

However, the component ^^^ = A{ri, rf) is perpendicular to the plane P and is therefore invariant under 5*0(2) 
rotations in that plane, so that it behaves like any other scalar function defined on P. It follows that the collection 
of integrals (27) is indeed the standard 2D Radon transform[2>2>\ of the scalar function A^^{-x.), x e P{y,r]), and 
hence can be inverted for the component A^^ using any numerical scheme for Radon inversion appropriate to the 
circumstances, such as filtered back-projection[34]. This produces the component A^^(x) for every point x in the 
plane P{y,ri). 

On repeating the same process for all planes parallel to P(y, J7) we reconstruct the component A^^ within the whole 
specimen. 

We now perform this procedure for the following six different choices of the vector rj: 

(28a) 



fll = 


= ei , ?72 = 


= e2 , 
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1723 = 
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/77 



mi = ^— • (28b) 

V'z yz \/2 

The scan-and-reconstruction cycle for the choices in (28a) produces components 

All , ^22 , A33 (29) 

of the anisotropy tensor. On the other hand, for the choices in (28b) we obtain the following result: Let us focus 
attention on the first orientation r/12, where the associated reconstruction gives us the tensor component ^(■'712, ''?12) 
everywhere within the object. But, due to the fact that the tensor is linear in its arguments, this component can be 
expressed in the form 

^(r?i2,r?i2) = -(^ii+^22)+^i2 , (30) 

where we have used the fact that Aij is symmetric and hence A21 = Ai2. Thus, having already reconstructed An 
and ^4.22 everywhere in the specimen, we can immediately compute A12 from the reconstructed values of A(rji2,rji2) 
using eq. (30). On repeating this process for the last two choices of t] in (28b) we see that all six tensor components 
of Aij can be reconstructed in this way. 

If the average dielectric constant e of the object is known we can compute the full dielectric tensor e^ immediately 
using the definition (8a). On the other hand, if e is not determined we can still use (8a) to write 

ey =e(Ajj +(5y) , (31) 

in other words, we can reconstruct eij up to a scale factor e; this may still produce valuable information about the 
internal structure of the dielectric material, see Fig. 1. 

VI. NUMERICAL EXAMPLE 

Here we present a numerical example of reconstruction of a single tensor component A^^ for a plane intersecting 
the object at an oblique angle. The polarization transformation data are obtained from an artificial stress model of a 
cylindrical bar with a circular cross-section which is subject to axial load and in turn bulges out in the middle section, 
see Fig. 1(a). 

Based on these artificial forward data we then employ our method and reconstruct the "normal" component A^^ on 
a plane passing through the center of the cylinder and making an angle of 22° with the symmetry axis, see Fig. 1(b). 
The original plot of A^^ in this plane is shown in Fig. 1(c); the reconstructed image in Fig. 1(d) has 254 x 254 pixels, 
assuming that 36 scans, one scan on every 5°, have been performed in the plane. It will be seen that the reconstruction 



contains artefacts which are typical of a Radon inversion; they can be reduced by increasing the number of scans, e.g., 
by performing one scan on each degree, for a total of 180 scans. The result in this case is almost indistinguishable 
from the original in Fig. 1(c) so that we have refrained from showing it. 

For the inverse Radon transformation, the Matlab function iradon has been used which utilizes a filtered back- 
projection algorithm [34]. 

VII. SUMMARY 

We have presented a novel way to reconstruct an arbitrarily inhomogeneous anisotropic dielectric tensor inside a 
transparent non-absorbing medium, under the conditions that the birefringence is a deviation from a homogeneous 
isotropic average, and that this deviation is weak. It was shown that the associated linearized inverse problem of 
reconstructing the dielectric tensor from polarization transformation data gathered by optical-tomographical means 
can be reduced to six scalar Radon inversions which allow the determination of the permittivity tensor completely. We 
also supplied a more rigorous derivation of the usual equations of integrated photoelasticity which define the inverse 
problem for the dielectric tensor; our exposition was based on a careful description of the various approximations that 
enter the derivation of the photoelasticity equations from Maxwell's equations in a material medium. 



Ackno\vledgements 

The authors acknowledge support from EPSRC grant GR/86300/01. 




(a) Cylindrical bar 
axially loaded 



(b) Oblique 
intersection 




(c) Original tensor 
component Ann 



(d) Reconstruction 

with 254 X 254 pixels 

and 36 scans 



FIG. 1: The reconstruction of the anisotropy tensor (8a) gives valuable information about the internal structure of the object, 
even if the average dielectric constant e, eq. (1), is not known. 
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